用 R 來玩空間分析：

本文目錄（點擊可跳轉）

[TOC]

---

library(magrittr)
library(data.table)

1. 基本地圖資料：套件 sp

套件 sp 提供了空間資料的 classes 和 methods，讓空間資料可以呈現 points 、 lines 、 polygons 和 grids 的特性

1-1. 空間資料的資料結構

2. 基本地圖資料處理：套件rgdal

Geospatial Data Abstraction Library(GDAL) 是由 Open Source Geospatial Fundation 所開發與維護的套件，而 rgdal 正是將其與 R 綁在一起的 R 語言套件。

2-1. 載入資料

setwd("~/Dropbox/practice/RGIS")
library(rgdal)

2-1-1. 載入經緯度資料文件並繪製地圖

TPE = read.table("data/TWN/tpedata.csv", header=TRUE, sep=",")
XCOOR = TPE$X
YCOOR = TPE$Y

#設定參考坐標系統
proj <- CRS("+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs")

利用 sp::SpatialPointsDataFrame 函式，將資料轉為可做進一步空間分析的 SpatialPointsDataFrame 資料格式

data.sp <- SpatialPointsDataFrame(coords=cbind(x=XCOOR,y=YCOOR), data=TPE, proj4string=proj)
# 看一下 data.sp 的 structure 長什麼樣子
data.sp %>% str

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  63 obs. of  3 variables:
##   .. ..$ WEEK: int [1:63] 4 42 31 8 50 33 13 26 24 23 ...
##   .. ..$ X   : int [1:63] 305197 306471 305461 305946 304515 305057 305413 303234 307774 307974 ...
##   .. ..$ Y   : int [1:63] 2763877 2764321 2764674 2764799 2765329 2765560 2765936 2766068 2766142 2766376 ...
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:63, 1:2] 305197 306471 305461 305946 304515 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ bbox       : num [1:2, 1:2] 296621 2763877 310825 2781440
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs"

畫成地圖來看看

plot(data.sp,pch=1, col="blue")

2-1-2. 載入 shp 地圖資料

地圖資料包含有數個檔案，以範例圖檔 country 為例，圖檔資料夾包括 country.dbf 、 country.prj 、 country.shp 和 country.shx，因此在讀入檔案時要先指向資料路徑，在 layer 參數中寫入圖檔的名字

World <- readOGR(dsn = "data/World", layer = "country")

## OGR data source with driver: ESRI Shapefile 
## Source: "data/World", layer: "country"
## with 251 features
## It has 11 fields

WorldCity <- readOGR(dsn = "data/World", layer = "CITIES")

## OGR data source with driver: ESRI Shapefile 
## Source: "data/World", layer: "CITIES"
## with 606 features
## It has 4 fields

# 看一下他們的 structure 長什麼樣子，因為結構龐大，所以看到第二層就好
World %>% str(max.level=2)

## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  251 obs. of  11 variables:
##   ..@ polygons   :List of 251
##   .. .. [list output truncated]
##   ..@ plotOrder  : int [1:251] 15 192 36 233 42 30 12 87 127 106 ...
##   ..@ bbox       : num [1:2, 1:2] -180 -90 180 83.6
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

WorldCity %>% str(max.level=2)

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  606 obs. of  4 variables:
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:606, 1:2] 33.1 40.6 30.5 150.8 56.2 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ bbox       : num [1:2, 1:2] -165.3 -53.2 177.1 78.2
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

可以看到 World 的資料格式是 SpatialPolygonsDataFrame 而 WorldCity 則是 SpatialPointsDataFrame

#看一下 圖檔的資料概況
summary(World)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##    min      max
## x -180 180.0000
## y  -90  83.6236
## Is projected: FALSE 
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
##    FIPS_CNTRY    GMI_CNTRY            CNTRY_NAME           SOVEREIGN  
##  AA     :  1   UMI    :  6   Afghanistan   :  1   United Kingdom: 16  
##  AC     :  1   ISR    :  3   Albania       :  1   France        : 13  
##  AF     :  1   REU    :  3   Algeria       :  1   United States : 13  
##  AG     :  1   NOR    :  2   American Samoa:  1   Australia     :  5  
##  AJ     :  1   SJM    :  2   Andorra       :  1   New Zealand   :  4  
##  AL     :  1   ABW    :  1   Angola        :  1   Norway        :  4  
##  (Other):245   (Other):234   (Other)       :245   (Other)       :196  
##    POP_CNTRY            SQKM_CNTRY         SQMI_CNTRY          CURR_TYPE  
##  Min.   :    -99999   Min.   :       2   Min.   :      1   Dollar   : 29  
##  1st Qu.:    111033   1st Qu.:     717   1st Qu.:    277   Franc    : 16  
##  Median :   3084641   Median :   61937   Median :  23914   CFA Franc: 13  
##  Mean   :  22696879   Mean   :  583961   Mean   : 225467   Pound    :  9  
##  3rd Qu.:  11316975   3rd Qu.:  350769   3rd Qu.: 135432   EC Dollar:  8  
##  Max.   :1281008318   Max.   :16851940   Max.   :6506534   (Other)  :151  
##                                                            NA's     : 25  
##    CURR_CODE   LANDLOCKED   COLOR_MAP 
##  USD    : 12   N:208      5      :35  
##  XOF    : 11   Y: 43      4      :34  
##  FRF    :  8              1      :33  
##  XCD    :  8              8      :32  
##  AUD    :  4              7      :31  
##  (Other):165              3      :30  
##  NA's   : 43              (Other):56

summary(WorldCity)

## Object of class SpatialPointsDataFrame
## Coordinates:
##               min      max
## coords.x1 -165.27 177.1302
## coords.x2  -53.15  78.2000
## Is projected: FALSE 
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Number of points: 606
## Data attributes:
##         NAME        COUNTRY      POPULATION       CAPITAL
##  Hyderabad:  2   US     : 49   Min.   :     -99   N:442  
##  San Jose :  2   Russia : 37   1st Qu.:  243744   Y:164  
##  Tripoli  :  2   China  : 36   Median :  710447          
##  Valencia :  2   Canada : 22   Mean   : 1411999          
##  Abidjan  :  1   India  : 20   3rd Qu.: 1536250          
##  Abu Zaby :  1   Brazil : 19   Max.   :23620000          
##  (Other)  :596   (Other):423

2-2. 參考坐標系統（Coordinate Reference Systems）

在網站 Spatial Reference 中，可以查到各種座標系統的數學公式，且為了各種地圖平台使用方便，包括了諸多不同格式，非常方便。

以台灣慣用的平面座標系統 TWD97-TM2 為例， sp 套件中處理座標系統的函式 CRS 必須使用 proj4 格式 點選之後可以看到 TWD97 座標系統的公式是：

+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +units=m +no_defs

• proj = tmerc 投影方法為橫麥卡托

• lon_0 = 121 以經度 121 為地圖中央線

• k = 0.9999 中央經線尺度放大倍數為0.9999 中央經線與圓柱面相切密合，所以尺度為 1 會造成圖面其它地方被放大。為讓尺度變化較為均勻，於是將投影座標乘以某一常數，讓中央經線的尺度略小於1，逐漸往兩側放大，到投影帶中間某一部分尺度約為 1，投影帶邊緣則略大於 1，這個乘常數，稱做「中央經線尺度」

• x_0 = 250000 橫座標平移量為250000 為避免讓中央經線西側座標出現負數，而將投影座標加一個常數，即為「橫座標平移量」。

• ellps = GRS80 橢球參數使用橢球體 GRS80 地球真實形狀其實不是正圓，也不是完美的橢圓，而是為赤道地帶略為膨脹、兩極地區略成扁平接近，因此參考橢球體即為最接近地球真實外型的規則幾何形狀。而 GRS80 是 1980 年國際大地測量學與地球物理學協會（International Union of Geodesy and Geophysics，簡稱為IUGG）公布之參考橢球體（GRS80），其橢球參數為： 長半徑 a = 6378137 公尺，扁率 f = 1/298.257222101

關於座標系統，是一個複雜的故事了，有空再來發一篇專文。

2-2-1. 查看並修改投影參數

我們可以修改參考坐標系統的內容。以下舉例：

World@proj4string 

## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

# Robinson projection
wrld.rob <- spTransform(World, CRS("+proj=robin")) 
plot(wrld.rob)

# conical equidistant projection
wrld.eqdc <- spTransform(World, CRS("+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=90 +lat_2=90 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"))
plot(wrld.eqdc)

# Azimuthal Equidistant projection
wrld.aeqd <- spTransform(World, CRS("+proj=aeqd +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"))
plot(wrld.aeqd)

3. 繪製面量圖：套件 GISTools

GISTools 提供許多協助繪製面量圖的工具

library(GISTools)

3-1. 檢視地圖資料的屬性表格

# @data 可以看到地圖的屬性表格
World@data %>% head

##   FIPS_CNTRY GMI_CNTRY          CNTRY_NAME           SOVEREIGN POP_CNTRY
## 0         AA       ABW               Aruba         Netherlands     67074
## 1         AC       ATG Antigua and Barbuda Antigua and Barbuda     65212
## 2         AF       AFG         Afghanistan         Afghanistan  17250390
## 3         AG       DZA             Algeria             Algeria  27459230
## 4         AJ       AZE          Azerbaijan          Azerbaijan   5487866
## 5         AL       ALB             Albania             Albania   3416945
##    SQKM_CNTRY SQMI_CNTRY CURR_TYPE CURR_CODE LANDLOCKED COLOR_MAP
## 0     182.926     70.628    Florin       AWG          N         1
## 1     462.378    178.524 EC Dollar       XCD          N         2
## 2  641869.188 247825.703   Afghani       AFA          Y         3
## 3 2320972.000 896127.312     Dinar       DZD          N         3
## 4   85808.203  33130.551     Manat      <NA>          Y         4
## 5   28754.500  11102.110       Lek       ALL          N         6

# 現在要針對屬性資料中的人口密度 POP_CNTRY 欄位來做面量圖
wrld.pop <- World@data$POP_CNTRY
# 處理錯誤資料
wrld.pop[which(wrld.pop==0)]=NA
wrld.pop[which(wrld.pop==-99999)]=NA
wrld.pop %<>% na.omit

3-2. 資料處理

3-2-1. Select by Attributes

這邊應用 magrittr 和 data.table 套件做資料選取， 當然也可以用更傳統的方式來做選取

# 讀取台灣地圖資料
TWN <- readOGR(dsn = "data/TWN", layer = "Townships", encoding="big5")

## OGR data source with driver: ESRI Shapefile 
## Source: "data/TWN", layer: "Townships"
## with 368 features
## It has 5 fields

TWN@data %<>% as.data.table 
plot(TWN)
# 選取 6 個直轄市各區，面積大於 90 平方公里的區
sel1 = TWN@data %>% 
  .[,.I[County=="臺北市" | County=="新北市" | County=="桃園市" | County=="臺中市" | County=="臺南市" | County=="高雄市"] ]
sel2 = TWN@data %>% .[,.I[Area>90]]
plot(TWN[intersect(sel1, sel2),] ,col="red", add=TRUE)

TWN[intersect(sel1, sel2),]@data

##     UNI_ID     Town County Area Code
##  1:      6   新店區 新北市  120    6
##  2:      9   三峽區 新北市  192    9
##  3:     19   石碇區 新北市  144   19
##  4:     20   坪林區 新北市  171   20
##  5:     25   雙溪區 新北市  146   25
##  6:     26   貢寮區 新北市  100   26
##  7:     29   烏來區 新北市  321   29
##  8:     99   東勢區 臺中市  117   99
##  9:    115   霧峰區 臺中市   98  115
## 10:    116   太平區 臺中市  121  116
## 11:    118   和平區 臺中市 1038  118
## 12:    199   安南區 臺南市  107  199
## 13:    204   白河區 臺南市  126  204
## 14:    207   東山區 臺南市  125  207
## 15:    216   七股區 臺南市  110  216
## 16:    225   楠西區 臺南市  110  225
## 17:    226   南化區 臺南市  172  226
## 18:    254   田寮區 高雄市   93  254
## 19:    262   旗山區 高雄市   95  262
## 20:    263   美濃區 高雄市  120  263
## 21:    264   六龜區 高雄市  194  264
## 22:    265   甲仙區 高雄市  124  265
## 23:    266   杉林區 高雄市  104  266
## 24:    267   內門區 高雄市   96  267
## 25:    268   茂林區 高雄市  194  268
## 26:    269   桃源區 高雄市  929  269
## 27:    270 那瑪夏區 高雄市  253  270
##     UNI_ID     Town County Area Code

3-2-2. Select by Locations

Flood <- readOGR(dsn = "data/TWN", layer = "flood50", encoding="big5")

## OGR data source with driver: ESRI Shapefile 
## Source: "data/TWN", layer: "flood50"
## with 5103 features
## It has 5 fields

Schools <- readOGR(dsn = "data/TWN", layer = "tpecity_school", encoding="big5")

## OGR data source with driver: ESRI Shapefile 
## Source: "data/TWN", layer: "tpecity_school"
## with 198 features
## It has 3 fields

# 先確保兩圖層的投影參數一樣
Flood@proj4string = CRS("+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs")

Schools@proj4string = CRS("+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs")
# 投影參數相同才能進行選取
Schools.Flood = Schools[Flood,] 

plot(Flood, col='cyan')
plot(Schools, col='yellow', pch=21, add=TRUE) 
plot(Schools.Flood, col='red', pch=4, add=TRUE) 

# 有哪些學校位處淹水區域？
(Schools.Flood@data)

##     TEXTSTRING   ID STUDENTS
## 1     士林國中    1       95
## 2     士林國小    2       93
## 8     明德國中   33      113
## 9     文林國小   39      213
## 13    富安國小   48       63
## 26    士東國小  300       52
## 42    大湖國小  681       92
## 52    百齡國小  730      163
## 61    新湖國小  837       24
## 67    民生國中  849       35
## 72    中山國中  863      200
## 73    五常國中  864       91
## 74    五常國小  865       55
## 75    吉林國小  882       91
## 78    劍潭國小  903      204
## 80    蓬萊國小  913      119
## 81    雙蓮國小  915       24
## 82    日新國小  916      201
## 83    民權國中  917      186
## 84    大同國小  918       25
## 87    太平國小  924       91
## 89    大橋國小  926      161
## 91    華江國小  979       99
## 98    龍山國小 1051       87
## 100   仁愛國小 1059       56
## 101   仁愛國中 1060      133
## 106   長春國小 1085       48
## 119   永吉國中 1199      101
## 123   木柵國中 1247      118
## 124   國光劇校 1248       68
## 131   興隆國小 1276       63
## 134   建安國小 1287       43
## 138   忠孝國小 1306      152
## 149   雙園國小 1362      110
## 174   興福國中 1575       85
## 175   溪口國小 1576      168

# 總共會影響多少學生？
sum(Schools.Flood@data$STUDENTS)

## [1] 3722

其中繪圖基本參數設定可參考此連結

3-2-3. Join Attributes

有人口資料 .csv 文件檔案，想要根據地區代號「unicode」， 和台灣地圖屬性資料中的地區代號「UNI_ID」互相 join 在一起

# 讀入人口資料 csv 檔案
Popn<-read.table("data/TWN/Townships_PopStat.csv", header=TRUE, sep=",") %>% data.table

head(Popn)

##    unicode towncode Male Female Age25 Age34 Age44 Age54 Age64 Age.L65
## 1:       1      101 6878   5146   205  1345  2433  3381  2747    1913
## 2:       2      102 4735   3464   122   948  1681  2277  1831    1340
## 3:       3      103 5267   3881   162  1110  1899  2488  2063    1426
## 4:       4      104 3057   2576   171   797  1204  1384  1162     915
## 5:       5      105 4768   3618   148  1115  2087  2414  1736     886
## 6:       6      106 3760   2766   128   631  1382  1826  1383    1176
##     Popn
## 1: 12024
## 2:  8199
## 3:  9148
## 4:  5633
## 5:  8386
## 6:  6526

head(TWN@data)

##    UNI_ID   Town County Area Code
## 1:      1 板橋區 新北市   23    1
## 2:      2 三重區 新北市   16    2
## 3:      3 中和區 新北市   20    3
## 4:      4 永和區 新北市    6    4
## 5:      5 新莊區 新北市   20    5
## 6:      6 新店區 新北市  120    6

# 要Join 之前，須確保要 join 的欄位名稱相同、型態也相同
# 這邊透過 data.frame 來修改欄位名稱與轉換型態
# 但 dplyr::left_join 只能使用 data frame 資料格式，所以還要轉回來
TWN@data %<>% .[,.(UNI_ID=UNI_ID %>% as.character,
                  Town, County, Area, Code)] %>% as.data.frame
Popn %<>% .[,.(UNI_ID=unicode %>% as.character,
               towncode, Male, Female, Age25, Age44, Age54, 
               Age64, Age.L65, Popn)] %>% as.data.frame
TWN@data<- dplyr::left_join(TWN@data, Popn)

## Joining by: "UNI_ID"

TWN@data %>% head

##   UNI_ID   Town County Area Code towncode Male Female Age25 Age44 Age54
## 1      1 板橋區 新北市   23    1      101 6878   5146   205  2433  3381
## 2      2 三重區 新北市   16    2      102 4735   3464   122  1681  2277
## 3      3 中和區 新北市   20    3      103 5267   3881   162  1899  2488
## 4      4 永和區 新北市    6    4      104 3057   2576   171  1204  1384
## 5      5 新莊區 新北市   20    5      105 4768   3618   148  2087  2414
## 6      6 新店區 新北市  120    6      106 3760   2766   128  1382  1826
##   Age64 Age.L65  Popn
## 1  2747    1913 12024
## 2  1831    1340  8199
## 3  2063    1426  9148
## 4  1162     915  5633
## 5  1736     886  8386
## 6  1383    1176  6526

3-2-4. Spatial Join

計算各台灣鄉鎮（Polygon）中所擁有的景點（Point）數量

Tour <- readOGR(dsn = "data/TWN", layer = "TWN_Tour")

## OGR data source with driver: ESRI Shapefile 
## Source: "data/TWN", layer: "TWN_Tour"
## with 943 features
## It has 3 fields

Tour@proj4string = proj
# selecting pts over polygon
Tour <- Tour[TWN, ] 

# 計算總數
# 這邊使用的函式 aggregate 是由 sp 套件所提供，專門處理 spatil objects
# aggregate 函式中所填入的 FUN 可以是計算個數的 length，也可以填入 mean, sd, max, min....或者自訂 function。
Tour_agg<-aggregate(x = Tour["Mark_ID"], by = TWN, FUN = length)  
summary(Tour_agg@data$Mark_ID)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.000   1.000   2.000   4.293   4.000  59.000     153

TWN$Tour_Count <- Tour_agg$Mark_ID
head(TWN@data) 

##   UNI_ID   Town County Area Code towncode Male Female Age25 Age44 Age54
## 1      1 板橋區 新北市   23    1      101 6878   5146   205  2433  3381
## 2      2 三重區 新北市   16    2      102 4735   3464   122  1681  2277
## 3      3 中和區 新北市   20    3      103 5267   3881   162  1899  2488
## 4      4 永和區 新北市    6    4      104 3057   2576   171  1204  1384
## 5      5 新莊區 新北市   20    5      105 4768   3618   148  2087  2414
## 6      6 新店區 新北市  120    6      106 3760   2766   128  1382  1826
##   Age64 Age.L65  Popn Tour_Count
## 1  2747    1913 12024          2
## 2  1831    1340  8199         NA
## 3  2063    1426  9148         NA
## 4  1162     915  5633         NA
## 5  1736     886  8386         NA
## 6  1383    1176  6526          7

# 處理 NA 值
sel = TWN@data %>% data.table %>% .[,.I[is.na(Tour_Count)]] 
TWN@data[sel,]$Tour_Count = 0
TWN@data %>% head

##   UNI_ID   Town County Area Code towncode Male Female Age25 Age44 Age54
## 1      1 板橋區 新北市   23    1      101 6878   5146   205  2433  3381
## 2      2 三重區 新北市   16    2      102 4735   3464   122  1681  2277
## 3      3 中和區 新北市   20    3      103 5267   3881   162  1899  2488
## 4      4 永和區 新北市    6    4      104 3057   2576   171  1204  1384
## 5      5 新莊區 新北市   20    5      105 4768   3618   148  2087  2414
## 6      6 新店區 新北市  120    6      106 3760   2766   128  1382  1826
##   Age64 Age.L65  Popn Tour_Count
## 1  2747    1913 12024          2
## 2  1831    1340  8199          0
## 3  2063    1426  9148          0
## 4  1162     915  5633          0
## 5  1736     886  8386          0
## 6  1383    1176  6526          7

3-2-5. Export to SHP

使用 rgdal 套件提供的 writeOGR 函式將剛才做出來的圖層存成 ESRI Shapefile 格式。 dsn 是資料夾路徑， layer 則填入圖檔名稱

writeOGR(TWN, dsn="data/output", layer="New_TWN", driver="ESRI Shapefile")

3-3. 繪製面量圖

3-3-1. 繪製面量圖與新增 legend、指北針與比例尺

新增 legend 之前要先設定 legend 要放在地圖的哪個位置， 執行 locator(1) 指令後點選一下剛才匯出的地圖， 會在 Console 中顯示剛才滑鼠點選的圖片 x 、 y 位置， 使用這個位置就可以指定 legend 和 指北針的位置

PDen = TWN$Popn/TWN$Area 
# locator(1)
# GISTools 提供的 auto.shading 函式可以將數值切成指定的分位數
shades <- auto.shading(PDen, cutter = quantileCuts, 
             n = 5, cols = brewer.pal(5, "RdBu"))
# 看一下 shades 的 structure ，可以看到其中已定義了 colors
shades %>% str

## List of 2
##  $ breaks: Named num [1:4] 15 29 53 99
##   ..- attr(*, "names")= chr [1:4] "20%" "40%" "60%" "80%"
##  $ cols  : chr [1:5] "#CA0020" "#F4A582" "#F7F7F7" "#92C5DE" ...
##  - attr(*, "class")= chr "shading"

# 使用 choropleth 函式，讀入剛才設定的分層 shades 就可以畫出地圖
choropleth(TWN, PDen, main="Taiwan Population Density")

choro.legend(343556.7,2590112, shades,cex=0.6, 
             fmt="%4.1f",title='Population Density (/km2)')
north.arrow(357257.5, 2644915, "N", len=8000, col="light gray")
map.scale(35288.13,2482789,100000,"100 km", 1, subdiv=1, tcol='black',scol='gray', sfcol='black')

3-3-2. 資料分級方法

在地圖學中，面量圖的資料如何分級是一門深奧的學問， 如果沒有根據資料的狀況做適當的分組， 在地圖著色之後可能會因為視覺效果影響，造成閱讀者判斷錯誤。

基本上有以下幾種分組方法： Defined Interval, Equal Interval, Quantile, Standard Deviation 和 Natural Break (Jenks)

我們可以透過 classInt 套件所提供的工具 classIntervals 來取得這些分級的結果

3-3-2-1. Defined Interval

自定義分級。例如 0 至 100 的數據， 定義 10 為間隔，那麼斷點則為 10、20、……、100， 分類數由間隔大小決定

library(classInt)
interval.DI = PDen %>% classIntervals(., 4, style="fixed",
                        fixedBreaks=c(0, 1000, 1500,2000 ,2557))

q <- cut(PDen, breaks= interval.DI$brks, include.lowest=T)
my_color = brewer.pal(9, "YlOrRd")[c(1,3,6,9)]

clr <- as.character(factor(q, labels = my_color))
plot(TWN, col = clr)

lenTEXT = 
c(paste0("<", interval.DI$brks[2]),    paste0(interval.DI$brks[2],"-",interval.DI$brks[3]),
paste0(interval.DI$brks[3],"-",interval.DI$brks[4]),  
paste0(interval.DI$brks[4],"-",interval.DI$brks[5]))

legend(legend = paste0(lenTEXT), fill = my_color, title="Pop Density", "topright")

3-3-2-2. Equal Interval

等距分級。固定組數，每組跨距相同。 例如 0～100 的數據，指定其分為 4 級，則間隔為 25

interval.EI = PDen %>% classIntervals(., 4, style="equal")

q <- cut(PDen, breaks= interval.EI$brks, include.lowest=T)
clr <- as.character(factor(q, labels = my_color))
plot(TWN, col = clr)

interval.EI$brks %<>% round(digits = 1)
lenTEXT = 
c(paste0("<", interval.EI$brks[2]),    paste0(interval.EI$brks[2],"-",interval.EI$brks[3]),
paste0(interval.EI$brks[3],"-",interval.EI$brks[4]),  
paste0(interval.EI$brks[4],"-",interval.EI$brks[5]))

legend(legend = paste0(lenTEXT), fill = my_color, title="Pop Density", "topright")

3-3-2-3. Quantile

Quantile 也就是前述 auto.shading 所提供的方法。 Quantile 分級適合用於線性分佈的資料。 但它不考慮數值大小，很可能將兩個大小相近的值分到不同的類別中， 也可能一樣的數據被分在不同的級距中。 在 Quantile 分級方法中，每一級距中的數量是一樣的

# 這邊也一樣用 classIntervals 示範一次 Quantile 分級方法
interval.QU = PDen %>% classIntervals(., 4, style="quantile")

q <- cut(PDen, breaks= interval.QU$brks, include.lowest=T)
clr <- as.character(factor(q, labels = my_color))
plot(TWN, col = clr)

interval.QU$brks %<>% round(digits = 1)
lenTEXT = 
c(paste0("<", interval.QU$brks[2]),    paste0(interval.QU$brks[2],"-",interval.QU$brks[3]),
paste0(interval.QU$brks[3],"-",interval.QU$brks[4]),  
paste0(interval.QU$brks[4],"-",interval.QU$brks[5]))

legend(legend = paste0(lenTEXT), fill = my_color, title="Pop Density", "topright")

3-3-2-4. Standard Deviation

標準差分級

interval.SD = PDen %>% classIntervals(., n=3, style="sd")

q <- cut(PDen, breaks= interval.SD$brks, include.lowest=T)
clr <- as.character(factor(q, labels = my_color[1:3]))
plot(TWN, col = clr)

interval.SD$brks %<>% round(digits = 1)
lenTEXT = 
c(paste0("<", interval.SD$brks[2]),    paste0(interval.SD$brks[2],"-",interval.SD$brks[3]),
paste0(interval.SD$brks[3],"-",interval.SD$brks[4])
)

legend(legend = paste0(lenTEXT), fill = my_color, title="Pop Density", "topright")

classIntervals 所提供的 standard deviation 分組方法，是設定好組數之後，系統自動計算要 mean±幾個sd 才能符合指定的組數。以上面設定的 3 組來說，各個斷點其實是 ( μ-5*σ-μ )、( μ-μ+5*σ)、( μ+5*σ- μ+10*σ) 這個以 5σ 為倍數的做法，只是為了達到指定組數為 3 的目的而設計， 因此如果要讓組數是傳統比較常看到的 μ±1*σ 或 μ±0.5*σ 可以用我寫的以下的函式 getSD：

getSD = function(data, nsd=1 ){
  
  u = mean(data)
  SD = sd(data)
  
  a = c(seq(from = 0, to = 20, by=nsd))
  
  stdUP = vector()
  for(i in 1:length(a)){
    stdUP[i]= u +SD*a[i+1]
    if( stdUP[i] > max(data) )
      break;
  }
  
  if (stdUP[length(stdUP)]>max(data)){  
    stdUP = stdUP[1:(length(stdUP)-1)]
  }else{
    stdUP=stdUP
  }
  
  n.upper = 1:length(stdUP)

  stdDW = vector()
  for(i in 1:length(a)){
    stdDW[i]= u -SD*a[i+1]
    if( stdDW[i] < min(data) )
      break;
  }
  
  if (stdDW[length(stdDW)] < min(data)){  
    stdDW = sort(stdDW[1:(length(stdDW)-1)])
  }else{
    stdDW=sort(stdDW)
  }
  
  n.lower = 1:length(stdDW)

  std=c(stdDW, u, stdUP)

  brks = paste0(c(-n.lower*nsd, 0,n.upper*nsd), "σ")
  brks[brks=="0σ"]="μ"
  return(list("std"=std, "brks"=brks))
}  

在getSD 函式中，nsd參數 即指定 σ 以 nsd 為倍數和 μ 做加減， 預設 nsd=1。 但要注意的是，當 nsd 設定得太小，可能會造成組數過多，不利閱讀。 在地圖學中建議組數不要超過 7 組。

std = getSD(PDen, nsd=1.5)$std
brks = getSD(PDen, nsd=1.5)$brks
q <- cut(PDen, breaks= std, include.lowest=T)
clr <- as.character(factor(q, 
                labels = brewer.pal(n = length(std)-1, "YlOrRd")))
plot(TWN, col = clr)

lenTEXT = 
c(paste0(brks[1], " - ", brks[2]),    
  paste0(brks[2], " - ", brks[3]),    
  paste0(brks[3], " - ", brks[4]),    
  paste0(brks[4], " - ", brks[5]),
  paste0(brks[5], " - ", brks[6]),
  paste0(brks[6], " - ", brks[7]),
  paste0(brks[7], " - ", brks[8]))


legend(legend = paste0(lenTEXT), 
       fill = brewer.pal(n = length(std)-1, "YlOrRd"), 
       title="Pop Density", "topright")

3-3-2-5. Natural Breaks (Jenks)

interval.NB = PDen %>% classIntervals(., 4, style="jenks")

q <- cut(PDen, breaks= interval.NB$brks, include.lowest=T)
clr <- as.character(factor(q, labels = my_color))
plot(TWN, col = clr)

interval.NB$brks %<>% round(digits = 1)
lenTEXT = 
c(paste0("<", interval.NB$brks[2]),    paste0(interval.NB$brks[2],"-",interval.NB$brks[3]),
paste0(interval.NB$brks[3],"-",interval.NB$brks[4]),  
paste0(interval.NB$brks[4],"-",interval.NB$brks[5]))

legend(legend = paste0(lenTEXT), fill = my_color, title="Pop Density", "topright")

3-3-3. 如何選擇適合的資料分級方法

如何判斷適合該資料的分級方法呢？ 應該先判斷該資料的分佈狀態。

hist(PDen)

boxplot(PDen)

可以從 histrogram 資料分佈中看到， PDen 屬於右偏態，也就是大多資料集中在左側 我們也可以從 boxplot 中看到，不管是下四分位數、中位數及上四分位數， 都非常集中在數值較低的區域

再來看看前述介紹幾種分級的方法的分組狀況：

3-3-3-1. Equal Interval

plot(interval.EI, pal=my_color, main="Equal Interval", xlab="PDen", ylab="")

Equal Interval 適合用在例如叫為人所知的資料，例如成績 (1-100分)、台灣氣溫 (0-40℃)、百分比(0-100%)等等。 但如果資料的分佈呈現 M 型分佈時，會導致某一塊區間出現空白。

而以 PDen 台灣人口密度此一資料而言，Equal Interval 的分級方法對閱讀者而言，這些分組斷點是「不具有任何統計意義」的， 因此雖然可以很好理解，我就是把人口密度從最小 0.929 到最大 2557 平均切成四份，但是對於解釋分組依據時，並沒有任何統計上的意義。

3-3-3-2. Quantile

plot(interval.QU, pal=my_color, main="Quantile", xlab="PDen", ylab="")

Quantile 分級方法，是將數據排名第25%、50%、75%的數值設為分組斷點， 也就是說每組中所擁有的資料數量是一樣的。例如 PDen 資料：

interval.QU

## style: quantile
##  [0.9,19.2)   [19.2,38)   [38,80.7) [80.7,2557] 
##          92          91          93          92

可以看到每組中的資料數量都是 92 個。 但是，它卻忽略了各組內的差異，例如最後一組 (80.7-2557) ，明明人口密度 80 人/km2 和 2557 人/km2 是超級巨大的差異，可是 Quantile 卻把他們編到為同一組了。

3-3-3-3. Standard Deviation

plot(interval.SD, pal=my_color, main="Standard Deviation", xlab="PDen", ylab="")

Standard Deviation 適合用在偏向常態分布的資料，其斷點具有很漂亮的統計意義。

3-3-3-4. Natural Breaking

plot(interval.NB, pal=my_color, main="Natural Breaks", xlab="PDen", ylab="")

Natural Break 所應用的就是最重要的分類原理：「讓組內差異最小」。

所謂的「差異」其實就是透過 variance 來表示，簡單的說， Natural Break 方法即是計算各種分類方法的 variance 總和，而 variance 總和最小的那種分類方法，就是最佳分類結果。 此方法計算量極大，詳細算法不在此贅述。

在示範資料 PDen 人口密度中，可以看到 Natural Break 可以呈現比較漂亮的結果，不會太集中於切割較低的數值，在高數值的呈現上，也不會讓組內的差異過大，是比較好的呈現方法。